
######################
## This R script replicates the analysis from "Increasing Intergovernmental Coordination to Fight Crime: Evidence from Mexico"
## The script includes results from the main paper and appendix
## Each code chunk is titled to correspond to the Table or Figure number in the published paper


# Set working directory ----
setwd("C:/Users/maa3411/Dropbox/Projects/Mando Unico/Data/")


# Start log ----
#sink("psrm_replication_log.txt")
library(logr)
tmp <- file.path("C:/Users/maa3411/Dropbox/Projects/Mando Unico/replication.log")
lf <- log_open(tmp)
options("logr.autolog" = TRUE, "logr.traceback" = FALSE)

# Packages -----
library(gsynth)
library(panelView)
library(dplyr)
library(cowplot)
library(ggplot2)
library(stargazer)
library(psych)

# Set seed -----
set.seed(1234)

# Load data ----

## annual data
load("psrm_cartel_data.RData")
# remove units that get rid of reform
yr_policial_synth <- mu_policial_annual[!(mu_policial_annual$idunico %in% c("11027","11010","11016","11033")),] 
# remove Leon, an outlier
yr_policial_synth <- yr_policial_synth[yr_policial_synth$idunico_num!=11020,]


## monthly data
load("psrm_crime_data.RData")
mt_policial_synth <- mu_policial_monthly
# exclude units that implement then remove reform
mt_policial_synth <- mt_policial_synth[!(mt_policial_synth$idunico %in% c("11016","11010","11027","11033","11005")),] 




# Treatment status ----

## FIGURE 1 (and A2): Treated status ----- 

panelview(presence_strength ~ policial, data = yr_policial_synth[yr_policial_synth$idunico_num!=11020,], 
          index = c("idunico","Year"), pre.post = TRUE, 
          by.timing = TRUE, ylab = "Municipality ID") 

## FIGURE A3: Treated status ----- 

panelview(presence_strength ~ Policial, data = mt_policial_synth, 
          index = c("idunico","time"), pre.post = TRUE, 
          by.timing = TRUE, ylab = "Municipality ID") 


## TABLE A1: List of municipality codes with municipality names, treatment status ----

treatment_list <- mt_policial_synth[mt_policial_synth$idunico!=11020,] %>% group_by(idunico, NOM_MUN) %>% summarise(treated = max(Policial))
treatment_list <- treatment_list[order(-treatment_list$treated),]
colnames(treatment_list) <- c("Municipality ID","Municipality name", "Ever treated")
a1 <- stargazer(as.data.frame(treatment_list), summary = F, header = F, 
          title = "List of municipalities in sample.",
          label = "tab:mun-list",
          out = "tex", rownames = F)
log_print("Table A1")
log_print(a1)


# Descriptive statistics -----

## TABLE A2 (TOP HALF) -----
year_summary_stats <- describe(yr_policial_synth[yr_policial_synth$idunico_num!=11020,c("policial", "presence_strength","co_num","war","log_pop","log_econin","state_political_vulnerability","federal_political_vulnerability","state_federal_political_vulnerability","TT_PERSO")],
                               fast = T)[,c(2:6)]
rownames(year_summary_stats) <- c("MUP", "Cartel presence strength", "Number of cartels", "Cartel war", "Log population", "Log economically inactive pop.",
                                  "Governor from rival party","President from rival party","Governor and president from rival party","Individuals in local public security")
a2_tophalf <- stargazer(year_summary_stats, summary = F, header = F, title = "Summary statistics for municipality-year analysis.",
          label = "tab:descriptive-stats-annual",
          out = "tex")
log_print("Table A2 - Top Half")
log_print(a2_tophalf)

## TABLE A2 (BOTTOM HALF) -----
month_summary_stats <- describe(mt_policial_synth[mt_policial_synth$idunico!=11020,c("Policial", "theft_violent_rate", "theft_nonviolent_rate", "hom_all_rate", "hom_ym_rate", "n.cell","n.weak","n.strong","war","log_pop","log_econin","state_political_vulnerability","federal_political_vulnerability","state_federal_political_vulnerability","TT_PERSO")],
                                fast = T)[,c(2:6)]
rownames(month_summary_stats) <- c("MUP", "Violent theft rate", "Non-violent theft rate", "Homicide rate", "Cartel-related homicide rate",
                                   "Number of cartel cells", "Number of weak cartels", "Number of strong cartels", "Cartel war",
                                   "Log population", "Log economically inactive pop.",
                                   "Governor from rival party","President from rival party","Governor and president from rival party","Individuals in local public security")
a2_bottomhalf <- stargazer(month_summary_stats, summary = F, header = F, title = "Summary statistics for municipality-month analysis.",
          label = "tab:descriptive-stats-month",
          out = "tex")
log_print("Table A2 - Bottom Half")
log_print(a2_bottomhalf)

# Main analysis ----

## Annual data ----
# cartel presence strength
out_strength <- gsynth(presence_strength ~ policial + log_pop + log_econin + state_political_vulnerability + federal_political_vulnerability + state_federal_political_vulnerability + TT_PERSO,# + mun_election_year, 
                       data = yr_policial_synth[yr_policial_synth$idunico_num!=11020,], EM = TRUE, 
                       index = c("idunico","Year"), force = "two-way", 
                       CV = TRUE, r = c(0, 5), se = TRUE, cl = "idunico",
                       inference = "parametric", nboots = 1000, 
                       parallel = T)

# number of cartels
out_number <- gsynth(co_num ~ policial + log_pop + log_econin + state_political_vulnerability + federal_political_vulnerability + state_federal_political_vulnerability + TT_PERSO,# + mun_election_year, 
                     data = yr_policial_synth[yr_policial_synth$idunico_num!=11020,], EM = TRUE, 
                     index = c("idunico","Year"), force = "two-way", 
                     CV = TRUE, r = c(0, 5), se = TRUE, cl = "idunico",
                     inference = "parametric", nboots = 1000, 
                     parallel = T)

# cartel war
out_war <- gsynth(war ~ policial + log_pop + log_econin + state_political_vulnerability + federal_political_vulnerability + state_federal_political_vulnerability + TT_PERSO,# + mun_election_year, 
                  data = yr_policial_synth[yr_policial_synth$idunico_num!=11020,], EM = TRUE, 
                  index = c("idunico","Year"), force = "two-way", 
                  CV = TRUE, r = c(0, 5), se = TRUE, cl = "idunico",
                  inference = "parametric", nboots = 1000, 
                  parallel = T)


## Monthly data ----

#theft violent
theft.viol.out2 <- gsynth(theft_violent_rate ~ Policial + n.cell + n.weak + n.strong + war + log_econin  + state_political_vulnerability + federal_political_vulnerability + state_federal_political_vulnerability + TT_PERSO, 
                          data = mt_policial_synth[(mt_policial_synth$Year>2010 & mt_policial_synth$idunico!=11020),], EM = TRUE, 
                          index = c("idunico","time"), force = "two-way", 
                          CV = TRUE, r = c(0, 5), se = TRUE, cl = "idunico",
                          inference = "parametric", nboots = 1000, 
                          parallel = T) 

#theft nonviolent
theft.nonviol.out2 <- gsynth(theft_nonviolent_rate ~ Policial + n.cell + n.weak + n.strong + war + log_econin + state_political_vulnerability + federal_political_vulnerability + state_federal_political_vulnerability + TT_PERSO, 
                             data = mt_policial_synth[(mt_policial_synth$Year>2010 & mt_policial_synth$idunico!=11020),], EM = TRUE, 
                             index = c("idunico","time"), force = "two-way", 
                             CV = TRUE, r = c(0, 5), se = TRUE, cl = "idunico",
                             inference = "parametric", nboots = 1000, 
                             parallel = T)

## Homicides
#all homs
hom.out2 <- gsynth(hom_all_rate ~ Policial + n.cell + n.weak + n.strong + war + log_econin + state_political_vulnerability + federal_political_vulnerability + state_federal_political_vulnerability + TT_PERSO, 
                   data = mt_policial_synth[(mt_policial_synth$Year<2021 & mt_policial_synth$idunico!=11020),], EM = TRUE, 
                   index = c("idunico","time"), force = "two-way", 
                   CV = TRUE, r = c(0, 5), se = TRUE, cl = "idunico",
                   inference = "parametric", nboots = 1000, 
                   parallel = T)
#homicides of young men
hom.ym.out2 <- gsynth(hom_ym_rate ~ Policial + n.cell + n.weak + n.strong + war + log_econin + state_political_vulnerability + federal_political_vulnerability + state_federal_political_vulnerability + TT_PERSO, 
                      data = mt_policial_synth[(mt_policial_synth$Year<2021 & mt_policial_synth$idunico!=11020),], EM = TRUE, 
                      index = c("idunico","time"), force = "two-way", 
                      CV = TRUE, r = c(0, 5), se = TRUE, cl = "idunico",
                      inference = "parametric", nboots = 1000, 
                      parallel = T)



# RESULTS: TABLES ----

### TABLE 1: Average ATT cartels ----

# empty models for reg table
a <- rnorm(100)
b <- rnorm(100)
c <- rnorm(100)
m1a <- lm(a ~ runif(100))
m2a <- lm(b ~ runif(100))
m3a <- lm(c ~ runif(100))

# add Avg. ATT to models
m1a$coefficients[2] <- out_strength$est.avg[1]
m2a$coefficients[2] <- out_number$est.avg[1]
m3a$coefficients[2] <- out_war$est.avg[1]

# Table 1
t1 <- stargazer(m1a, m2a, m3a, header = F, omit = "Constant", 
          add.lines = list(c("Municipality FE","Yes","Yes","Yes"),
                           c("Year FE","Yes","Yes","Yes"),
                           c("Unobserved factors","1","1","1"),
                           c("Period","2000-2021","2000-2021","2000-2021","2000-2021"),
                           c("Observations","726","726","726"),
                           c("Treated Muns","10","10","10"),
                           c("Control Muns","23","23","23")),
          se = list(c(1,out_strength$est.avg[2]),c(1,out_number$est.avg[2]),c(1,out_war$est.avg[2])),
          p = list("Police Reform"=c(1,out_strength$est.avg[5],out_number$est.avg[5],out_war$est.avg[5])),
          omit.stat = c("all"),
          covariate.labels = c("Police Reform"),
          label = "tab:main-annual",
          dep.var.labels = c("Cartel strength","Cartel number","Cartel war"),
          out = "tex")
log_print("Table 1")
log_print(t1)


### TABLES A4, A5, and A6: ATT cartels -----

strength_att <- out_strength$est.att[15:22,]
number_att <- out_number$est.att[15:22,]
war_att <- out_war$est.att[15:22,]

a4 <- stargazer(strength_att, summary = F, header = F, 
          title = "ATT effect of increased intergovernmental coordination on cartel strength of presence per treatment period.",
          label = "tab:strength-att", digits = 3,
          out = "tex")
a5 <- stargazer(number_att, summary = F, header = F, 
          title = "ATT effect of increased intergovernmental coordination on number of cartels per treatment period.",
          label = "tab:number-att", digits = 3,
          out = "tex")
a6 <- stargazer(war_att, summary = F, header = F, 
          title = "ATT effect of increased intergovernmental coordination on cartel wars per treatment period.",
          label = "tab:war-att", digits = 3,
          out = "tex")
log_print("Table A4")
log_print(a4)
log_print("Table A5")
log_print(a5)
log_print("Table A6")
log_print(a6)


## TABLE 2: Average ATT - violence and crime ----

# empty models for reg table
d <- rnorm(100)

m1m2 <- lm(a ~ runif(100))
m2m2 <- lm(b ~ runif(100))
m3m2 <- lm(c ~ runif(100))
m4m2 <- lm(d ~ runif(100))

# add Avg. ATT to models
m1m2$coefficients[2] <- theft.viol.out2$est.avg[1]
m2m2$coefficients[2] <- theft.nonviol.out2$est.avg[1]
m3m2$coefficients[2] <- hom.out2$est.avg[1]
m4m2$coefficients[2] <- hom.ym.out2$est.avg[1]

## Table 2
t2 <- stargazer(m1m2, m2m2, m3m2, m4m2, header = F,
          add.lines = list(c("Municipality FE","Yes","Yes","Yes","Yes"),
                           c("Year FE","Yes","Yes","Yes","Yes"),
                           c("Unobserved factors","0","1","2","1"),
                           c("Period","1/2011-12/2021","1/2011-12/2021","1/2000-12/2020","1/2000-12/2020"),
                           c("Observations","4356","4356","8316","8316"),#
                           c("Treated Muns","11","11","10","10"),
                           c("Control Muns","22","22","23","23")),
          se = list(c(1,theft.viol.out2$est.avg[2]),c(1,theft.nonviol.out2$est.avg[2]),c(1,hom.out2$est.avg[2]),c(1,hom.ym.out2$est.avg[2])),
          p = list("Police Reform"=c(1,theft.viol.out2$est.avg[5],theft.nonviol.out2$est.avg[5],hom.out2$est.avg[5],hom.ym.out2$est.avg[5])),
          omit.stat = c("all"), omit = "Constant", 
          covariate.labels = c("Police Reform"),
          label = "tab:main-monthly-rate",
          dep.var.labels = c("Violent theft","Non-violent theft","All homicides","Homicides, young men"),
          out = "tex")
log_print("Table 2")
log_print(t2)


## TABLES A7, A8, A9, A10: ATT Crimes -----

vtheft_att1 <- theft.viol.out2$est.att[42:112,]
nvtheft_att1 <- theft.nonviol.out2$est.att[42:112,]
hom_att1 <- hom.out2$est.att[174:244,]
homym_att1 <- hom.ym.out2$est.att[174:244,]

vtheft_att2 <- theft.viol.out2$est.att[42:112,]
nvtheft_att2 <- theft.nonviol.out2$est.att[42:112,]
hom_att2 <- hom.out2$est.att[174:244,]
homym_att2 <- hom.ym.out2$est.att[174:244,]


a7 <- stargazer(vtheft_att2, summary = F, header = F, 
          title = "ATT effect of increased intergovernmental coordination on violent theft rates per treatment period.",
          label = "tab:vtheft-att", digits = 3,
          out = "tex")
a8 <- stargazer(nvtheft_att2, summary = F, header = F, 
          title = "ATT effect of increased intergovernmental coordination on non-violent theft rates per treatment period.",
          label = "tab:nvtheft-att", digits = 3,
          out = "tex")
a9 <- stargazer(hom_att2, summary = F, header = F, 
          title = "ATT effect of increased intergovernmental coordination on homicide rates per treatment period.",
          label = "tab:hom-att", digits = 3,
          out = "tex")
a10 <- stargazer(homym_att2, summary = F, header = F, 
          title = "ATT effect of increased intergovernmental coordination on cartel-related homicide rates per treatment period.",
          label = "tab:homym-att", digits = 3,
          out = "tex")

log_print("Table A7")
log_print(a7)
log_print("Table A8")
log_print(a8)
log_print("Table A9")
log_print(a9)
log_print("Table A10")
log_print(a10)



# RESULTS: FIGURES ----

## FIGURE 2: Cartel results ----

# strength
p1a <- plot(out_strength, type = "counterfactual", raw = "none", main= "", ylab = "Cartel strength", legendOff = T)
#plot(out, type = "counterfactual", raw = "all", main="Cartel strength of presence (0-3)")
p2a <- plot(out_strength, main= "")

# number of cartels
p3a <- plot(out_number, type = "counterfactual", raw = "none", main= "", ylab = "Number of cartels", legendOff = T)
p4a <- plot(out_number, main= "")

# (create legends)
plot_legend1 <- ggplot(yr_policial_synth[yr_policial_synth$idunico %in% c(11001,11002),], aes(x=Year,y=log_pop, group = idunico, color = idunico)) + 
  geom_line(lwd = 2) + theme(legend.position="bottom", legend.key = element_rect(fill = "white")) + scale_color_manual(values = c("black", "darkblue"),name="", breaks = c("11001","11002"), labels = c("Treated Average","Estimated Y(0) Average"))
plot_legend2 <- ggplot(yr_policial_synth[yr_policial_synth$idunico %in% c(11001,11002),], aes(x=Year,y=log_pop, group = idunico, color = idunico)) + 
  geom_line(lwd = 2) + theme(legend.position="bottom", legend.key = element_rect(fill = c("grey"))) + scale_color_manual(values = c("black", "grey"),name="", breaks = c("11001","11002"), labels = c("Estimated ATT","95% Confidence Interval"))

legend1 <- get_legend(plot_legend1)
legend2 <- get_legend(plot_legend2)

# cartel war
p5a <- plot(out_war, type = "counterfactual", raw = "none", main= "", ylab = "Cartel war", legendOff = T)
p6a <- plot(out_war, main= "")


## plot all together 
plot_grid(p1a, p2a, p3a, p4a, p5a, p6a, legend1, legend2,
          nrow = 4,
          labels = c("(A1)","(A2)","(B1)","(B2)","(C1)","(C2)","",""),
          label_size = 12,
          label_y = .12,
          rel_heights = c(1,1,1,.3)
)




## FIGURE 3: Violence and crime -----

# violent theft
p1m2 <- plot(theft.viol.out2, type = "counterfactual", raw = "none", main="", ylab = "Theft, violent\nper 100,000", legendOff = T)
p2m2 <- plot(theft.viol.out2, main= "")
p1m2b <- plot(theft.viol.out2, type = "counterfactual", raw = "none", main="", 
              ylab = "Theft, violent\nper 100,000", legendOff = T, xlim = c(-75,100))

# nonviolent theft
p3m2 <- plot(theft.nonviol.out2, type = "counterfactual", raw = "none", main="", ylab = "Theft, non-violent\nper 100,000", legendOff = T)
p4m2 <- plot(theft.nonviol.out2, main= "")
p3m2b <- plot(theft.nonviol.out2, type = "counterfactual", raw = "none", main="", 
              ylab = "Theft, non-violent\nper 100,000", legendOff = T, xlim = c(-75,100))

# all homicides
p7m2 <- plot(hom.out2, type = "counterfactual", raw = "none", main="", ylab = "Homicides, all\nper 100,000", legendOff = T)
p8m2 <- plot(hom.out2, main= "")
p7m2b <- plot(hom.out2, type = "counterfactual", raw = "none", main="", 
              ylab = "Homicides, all\nper 100,000", legendOff = T, xlim = c(-100,100))
p8m2b <- plot(hom.out2, main= "", xlim = c(-100,100))

# homicides, young men
p9m2 <- plot(hom.ym.out2, type = "counterfactual", raw = "none", main="", ylab = "Homicides, young men\nper 100,000", legendOff = T)
p10m2 <- plot(hom.ym.out2, main= "")
p9m2b <- plot(hom.ym.out2, type = "counterfactual", raw = "none", main="", 
              ylab = "Homicides, young men\nper 100,000", legendOff = T, xlim = c(-100,100))
p10m2b <- plot(hom.ym.out2, main= "", xlim = c(-100,100))

## plot all together 
plot_grid(p1m2b, p2m2, p3m2b, p4m2, p7m2b, p8m2b, p9m2b, p10m2b, legend1, legend2,
          nrow = 5,
          labels = c("(A1)","(A2)","(B1)","(B2)","(C1)","(C2)","(D1)","(D2)","",""),
          label_size = 12,
          label_y = .12,
          rel_heights = c(1,1,1,1,.3)
)





# Outcome trends ----

## FIGURE A5 -----

id_policial <- unique(yr_policial_synth$idunico[yr_policial_synth$policial==1])
yr_trends_data <- yr_policial_synth[yr_policial_synth$idunico_num!=11020,]
yr_trends_data$ever_policial <- NA
yr_trends_data$ever_policial[yr_trends_data$idunico %in% id_policial] <- 1
yr_trends_data$ever_policial[!(yr_trends_data$idunico %in% id_policial)] <- 0

annual_trends <- yr_trends_data %>% group_by(Year,ever_policial) %>% 
  summarise(strength = mean(presence_strength), co_num = mean(co_num), war = mean(war))
annual_trends$ever_policial <- as.factor(annual_trends$ever_policial)

annual_trends %>%
  ggplot( aes(x=Year, y=strength, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Cartel Strength") + 
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="none") + geom_vline(xintercept = 2014)
annual_trends %>%
  ggplot( aes(x=Year, y=co_num, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Cartel Number") + 
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="none") + geom_vline(xintercept = 2014)
annual_trends %>%
  ggplot( aes(x=Year, y=war, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Cartel Wars") + 
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="bottom") + geom_vline(xintercept = 2014)

## FIGURE A6 ----
mt_id_policial <- unique(mt_policial_synth$idunico[mt_policial_synth$Policial==1])
mt_trends_data <- mt_policial_synth[mt_policial_synth$idunico!=11020,]
mt_trends_data$ever_policial <- NA
mt_trends_data$ever_policial[mt_trends_data$idunico %in% mt_id_policial] <- 1
mt_trends_data$ever_policial[!(mt_trends_data$idunico %in% mt_id_policial)] <- 0

annual_trends_m <- mt_trends_data %>% group_by(time,ever_policial) %>% 
  summarise(vtheft = mean(theft_violent,na.rm = T), nvtheft = mean(theft_nonviolent,na.rm = T),
            hom = mean(hom_all,na.rm = T), homym = mean(hom_ym,na.rm = T))
annual_trends_m$ever_policial <- as.factor(annual_trends_m$ever_policial)

annual_trends_m %>%
  ggplot( aes(x=time, y=vtheft, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Violent Theft") + xlab("Time") +
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="none") + 
  geom_vline(xintercept = 175) + xlim(130,270)
annual_trends_m %>%
  ggplot( aes(x=time, y=nvtheft, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Non-Violent Theft") + xlab("Time") +
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="none") + 
  geom_vline(xintercept = 175) + xlim(130,270)
annual_trends_m %>%
  ggplot( aes(x=time, y=hom, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Homicides") + xlab("Time") +
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="none") + 
  geom_vline(xintercept = 175) + xlim(1,270) 
annual_trends_m %>%
  ggplot( aes(x=time, y=homym, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Homicides, Young Men") + xlab("Time") +
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="none") + 
  geom_vline(xintercept = 175) + xlim(1,270) 

# get legend
library(ggpubr)
level_legend <- annual_trends_m %>%
  ggplot( aes(x=time, y=homym, group=ever_policial, color=ever_policial)) +
  geom_line(aes(linetype=ever_policial), size = 1.5) + theme_bw() +
  theme_bw() + ylab("Homicides, yYoung Men") + xlab("Time") +
  scale_color_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  scale_linetype_discrete(labels = c("No", "Yes"), breaks = c(0, 1), name = "Implemented reform") +
  theme(text = element_text(size = 20),legend.position="bottom") + 
  geom_vline(xintercept = 175) + xlim(1,270) 
as_ggplot(get_legend(level_legend))



# End log ----
#sink()
log_close()
writeLines(readLines(lf))
